clear
%% Problem
%   'EqualMinima','Himmelblau','Sixhump','ModifiedRastrigin','Vincent',
%   'Shubert','Composition','IncreasingMinima','Rastrigin','Schaffer'.
ProblemSet.FuncType = 'Vincent';
ProblemSet.dimension = 2;
% ProblemSet.k = [3,3,3];
% ProblemSet.BaseFunc = 3;

[d, x_bound, ObjFunc, OptSet, precision_init, precision, ~, MaximalBudget, FileName]...
    = problem_setting(ProblemSet);

MaximalBudget = 10000;
iepsilon = 4;
epsilon = [0.1,0.01,0.001,0.0001];
OptSet.epsilon = epsilon(iepsilon);

%% problem
Problem.Dimension = d;
Problem.Domain = reshape(x_bound,[1,d*2]);
Problem.Sampling = @(n, region)Sampling_BoxConstraint(n, region, ObjFunc);
NewRegionNum = 2;
MaxPartitionDepth_d = ceil(log((x_bound(2,:)-x_bound(1,:))./precision_init)./log(NewRegionNum));
MaxPartitionDepth = sum(MaxPartitionDepth_d);
Problem.Partition = @(Region,SampleSet)Partition_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
%% Algorithm
AlgorithmM.n0 = 6;
AlgorithmM.QuantileLevel = 0.3;
AlgorithmM.NewBudget = 3;
AlgorithmM.MaxSampleSize = 15;%5*(d+1);
AlgorithmM.StopCriteria = [1,MaximalBudget];
ClearRadius0 = (x_bound(2,:)-x_bound(1,:))./(NewRegionNum.^(MaxPartitionDepth_d));
ClearRadius = 2*min(ClearRadius0); % the maximum distance within a non-partitionable region
AlgorithmM.radius = ClearRadius;
AlgorithmM.local_search = @(starts,step,MaximalBudget)LocalSearch_test(ObjFunc,x_bound,starts,step,precision(iepsilon)/2,OptSet,MaximalBudget);
%% OPTIMIZATION
rep = 100;
OptObtained = zeros(rep,1);
TotalBudget1 = zeros(rep,1);
TotalBudget2 = zeros(rep,1);
OptNoFound = zeros(rep,1);
Times = zeros(rep,1);
OptRemain = cell(rep,1);
%%
for i=1:rep
    disp(['epsilon: ',num2str(epsilon(iepsilon)),', replication: ',num2str(i),'.'])
    tic
    [optima, ~, SampleSet, ~, ls_budget, OptSetRemain]...
        = prsmmo_ls_test( Problem, AlgorithmM, OptSet );
    Times(i,1) = toc;
    disp(['time used: ',num2str(Times(i,1)),'.'])
    OptObtained(i,1) = size(optima,1);
    TotalBudget2(i,1) = ls_budget;
    TotalBudget1(i,1) = size(SampleSet,1) - ls_budget;
    OptRemain{i,1} = OptSetRemain;
    OptNoFound(i,1) = size(OptSetRemain,1);
end
%%
clear i iepsilon OptSetRemain SampleSet ls_budget optima
FileName = ['Vincent2D_a',num2str(AlgorithmM.QuantileLevel),'_N',num2str(MaximalBudget),'.mat'];
save(FileName)
a=[mean(TotalBudget1)',mean(TotalBudget2)',mean(TotalBudget1+TotalBudget2)',var(TotalBudget1+TotalBudget2)'];
%% figures
range = [0,0.6,0.8,1];
OptNum = size(OptSet.sol,1);
FoundPerc = ones(OptNum,1).*rep;
for i = 1:rep
    id = ismember(OptSet.sol,OptRemain{i,1},'rows');
    FoundPerc(id) = FoundPerc(id) - 1;
end
FoundPerc = FoundPerc./rep;
FoundPerc_view = reshape(FoundPerc,[sqrt(OptNum),sqrt(OptNum)])

figure()
hold on
rectangle('Position',[Problem.Domain([1,3]),Problem.Domain([2,4])-Problem.Domain([1,3])])
id = find((FoundPerc>=range(3)));
plot(OptSet.sol(id,1),OptSet.sol(id,2),'o','MarkerEdgeColor',[0.8500 0.3250 0.0980],'MarkerSize',5,'LineWidth',1.5)
id = find((FoundPerc<range(3)).*(FoundPerc>=range(2)));
plot(OptSet.sol(id,1),OptSet.sol(id,2),'^','MarkerEdgeColor',[0.9290 0.6940 0.1250],'MarkerSize',5,'LineWidth',1.5)
id = find((FoundPerc<range(2)).*(FoundPerc>=range(1)));
plot(OptSet.sol(id,1),OptSet.sol(id,2),'+','MarkerEdgeColor',[0 0.4470 0.7410],'MarkerSize',6,'LineWidth',1.5)
hold off
title('','Interpreter','latex','String',['$\alpha = ',num2str(AlgorithmM.QuantileLevel),'$'])
xlabel('','Interpreter','latex','String','$x_1$')
ylabel('','Interpreter','latex','String','$x_2$')
legend([num2str(range(3)),' \leq Pr \leq ',num2str(range(4))],[num2str(range(2)),' \leq Pr < ',num2str(range(3))],[num2str(range(1)),' \leq Pr < ',num2str(range(2))])
xlim(x_bound(:,1)')
ylim(x_bound(:,2)')
set(gca,'Fontname','Times New Roman','FontSize',16);
